close all; clear all; clc; 


psi1=0.3341;psi2=1.78;

psi1b=0.05;psi2b=0.02;

a = load('agrid2.txt');
a = exp(a);

%plot new graph


if (a(end)>psi1b/(2*psi2b));
    display(['problem' psi1b/(2*psi2b) a(end) ]);
    dbstop;
end;

udef = @(a) max(0,psi1+psi2*log(a));

% choose x
x1=0.9;
x2 = 1.1;
x3=0.7;

y1=udef(x1);
y2=udef(x2);
y3=udef(x3)-0.21;

%options = optimoptions('fsolve','Display','iter','Algorithm','trust-region');
%udef_param = @(cc) [cc(1)*x1+cc(2).*x1.^2-y1;cc(1)*x2+cc(2).*x2.^2-y2;cc(1)*x3+cc(2).*x3.^2-y3];

%[yy,fval,exitflag]=fsolve(udef_param,[5,1],options);


yy(1)=0.33;
yy(2)=1.38;
psi1b=yy(1);
psi2b=yy(2);

%exitflag

udefb= @(a) max(0,psi1b.*a+psi2b.*log(a));

figure;
plot(a,udef(a),'r',a,udefb(a),'b');
